##
## Replication data for JLC article
##

library(simcf)
library(ggplot2)
library(memisc)
#load mtable command after memisc
options(scipen=120,digits=2)

##
## Replicating Table 2
##

fates1 <- read.csv("replication_data_1.csv")

# Separate data sets are not required for modeling, but are useful in later stages like getting subset mean values for covariates and doing cross validation

form1 <- cname ~ year.x + ptfate2 + ptfateD + SLJI 
data1 <- extractdata(form1,fates1,na.rm=T)
form2 <- cname ~ year.x + ptfate2 + ptfateD + SLJI + chga_demo + Slog_gdp			
data2 <- extractdata(form2,fates1,na.rm=T)
form3 <- cname ~ year.x + ptfate2 + ptfateD + SLJI + chga_demo + Slog_gdp +
 	major_war + Syears_democ + entry + dpi_system	
data3 <- extractdata(form3,fates1,na.rm=T) #model 4 uses the same data as model 3


# TABLE 2 (logit results):

res.model1 <- glm(ptfateD ~ SLJI,data1,family=binomial)
res.model2 <- glm(ptfateD ~ SLJI + chga_demo + Slog_gdp,data2,family=binomial)
res.model3 <- glm(ptfateD ~ SLJI + chga_demo + Slog_gdp  + major_war +
			Syears_democ + as.factor(entry) + dpi_system,
			data3,family=binomial)		
res.model4 <- glm(ptfateD ~ SLJI + chga_demo,data3,family=binomial)

m1 <- mtable(res.model1,res.model2,res.model3,res.model4,			
		summary.stats=c("BIC","N"))			
m2 <- (m1 <- relabel(m1,
		"Ranef - cname" = "Random effect: country",
		"SLJI" = "Judicial independence",
		"chga_demo: 1. Democracy/0. Dictatorship" = "Democracy",
		"Slog_gdp" = "log(GDP/capita)",
		"major_war" = "Conflict",
		"Swdi_gdpgr" = "Economic growth",
		"Syears_democ" = "Years democratic",
		"as.factor(entry): 1/0" = "Entry by irregular means",
		"as.factor(entry): 2/0" = "Entry by foreign imposition",
		"dpi_system: 1. Strong president elected by assembly/0. Direct presidential" =
		"President elected by assembly",
		"dpi_system: 2. Parliamentary/0. Direct presidential" = "Parliamentary system"
		))
m2


##
## Replicating Figure 3
##

##
# MODEL ONE
##

myModel <- res.model1

nSims <- 1000

someScenarios <- expand.grid(1,
	seq(-0.4,0.6,0.01))

set.seed(1)	
simDraws <- mvrnorm(nSims,coef(myModel),vcov(myModel))
simYhats <- simDraws %*% t(someScenarios)
simYhats <- plogis(simYhats)

dim(simYhats)

predictionFrame <- data.frame(someScenarios,t(simYhats))

bulk <- apply(predictionFrame[1:101,3:1002],1,sort) 
# sort the predictions by their values, allowing us to remove the tails
bulk <- t(bulk) 
# needed to transpose the matrix
bulk2 <- bulk[,27:975]
# remove the tails
bulk3 <- cbind(predictionFrame[,1:2],bulk2)
# bind with the variable values from the prediction

library(reshape2)
longFrame <- melt(predictionFrame, id.vars = colnames(someScenarios))
head(longFrame)

# Below only if you want the mean of the simulations plotted as a line:
test1 <- seq(-0.4,0.6,0.01)
test2 <- tapply(longFrame$value,longFrame$Var2,mean)
variable <- rep(1,101)
test <- data.frame(test1,test2,variable)

library(ggplot2)
zp1 <- ggplot (data = longFrame,
		aes(x = Var2, y = value, group= paste(variable)))		
zp1 <- zp1 + geom_line(alpha=I(1/sqrt(nSims)))				
zp1 <- zp1 + geom_line(x=test1,y=test2,data=test,
		colour="white",size=1) # only for line at mean
zp1 <- zp1 + scale_x_continuous('',
		breaks=c(-0.4,-0.2,0,0.2,0.4,0.6),labels=c(
		"0.0","0.2","0.4","0.6","0.8","1.0")) #because "SLJI" is mean-centered
#zp1 <- zp1 + xlim(0,1)
zp1 <- zp1 + scale_y_continuous('Pr(Unpunished)',limits=c(0.2,1))
#zp1 <- zp1 + theme(plot.title="")
zp1 <- zp1 + theme_bw()
zp1 <- zp1 + theme(plot.title=element_text(size=25))
zp1 <- zp1 + theme(axis.title.y = element_text(size=24,angle=90),
          axis.text.y=element_text(size=20))
zp1 <- zp1 + theme(axis.title.x = element_text(size=18,vjust=0),
          axis.text.x  = element_text(size=20))
zp1          

##
# MODEL TWO
##
res.model2 <- glm(ptfateD ~ SLJI + chga_demo + Slog_gdp,data2,family=binomial)
democ <- mean(as.numeric(data2$chga_demo))-1 #-1 because as.numeric() changes it to 1-2 instead of 0-1, which is what is done in the model estimation
gdp <- mean(data2$Slog_gdp)

myModel <- res.model2

nSims <- 1000

someScenarios <- expand.grid(1,
	seq(-0.4,0.6,0.01),
	democ,gdp)
set.seed(1)
simDraws <- mvrnorm(nSims,coef(myModel),vcov(myModel))
simYhats <- simDraws %*% t(someScenarios)
simYhats <- plogis(simYhats)

dim(simYhats)

predictionFrame <- data.frame(someScenarios,t(simYhats))

# quantile(predictionFrame[1,3:1002],probs=c(0.025,0.975))
# determine the values empirically (I wanted to double check
# to make sure just doing the bulk of the dist was fine)

bulk <- apply(predictionFrame[1:101,5:1004],1,sort) 
# sort the predictions by their values, allowing us to remove the tails
bulk <- t(bulk) 
# needed to transpose the matrix
bulk2 <- bulk[,27:975]
# remove the tails
bulk3 <- cbind(predictionFrame[,1:4],bulk2)
# bind with the variable values from the prediction

longFrame <- melt(predictionFrame, id.vars = 
	colnames(someScenarios))
head(longFrame)

# Below only if you want the mean of the simulations plotted as a line:
test1 <- seq(-0.4,0.6,0.01)
test2 <- tapply(longFrame$value,longFrame$Var2,mean)
variable <- rep(1,101)
test <- data.frame(test1,test2,variable)

zp1 <- ggplot (data = longFrame,
		aes(x = Var2, y = value, group= paste(variable)))		
zp1 <- zp1 + geom_line(alpha=I(1/sqrt(nSims)))				
zp1 <- zp1 + geom_line(x=test1,y=test2,data=test,
		colour="white",size=1) # only for line at mean
zp1 <- zp1 + scale_x_continuous('Judicial independence',
		breaks=c(-0.4,-0.2,0,0.2,0.4,0.6),labels=c(
		"0.0","0.2","0.4","0.6","0.8","1.0"))
#zp1 <- zp1 + xlim(0,1)
zp1 <- zp1 + scale_y_continuous('',limits=c(0.2,1))
#zp1 <- zp1 + theme(plot.title="")
zp1 <- zp1 + theme_bw()
zp1 <- zp1 + theme(plot.title=element_text(size=25))
zp1 <- zp1 + theme(axis.title.y = element_text(size=24,angle=90),
          axis.text.y=element_text(size=20))
zp1 <- zp1 + theme(axis.title.x = element_text(size=24,vjust=0),
          axis.text.x  = element_text(size=20))
zp1          


##
# MODEL THREE
##
res.model3 <- glm(ptfateD ~ SLJI + chga_demo + Slog_gdp + major_war + 
		Syears_democ + as.factor(entry) + dpi_system,
		data3,family=binomial)

democ <- mean(as.numeric(data3$chga_demo))-1
gdp <- mean(data3$Slog_gdp)
conflict <- mean(data3$major_war)
years <- mean(data3$Syears_democ)

myModel <- res.model3
nSims <- 1000
someScenarios <- expand.grid(1,
	seq(-0.4,0.6,0.01),
	democ,gdp,conflict,years,0,0,0,0) # entry and executive set to reference categories
set.seed(1)	
simDraws <- mvrnorm(nSims,coef(myModel),vcov(myModel))
simYhats <- simDraws %*% t(someScenarios)
simYhats <- plogis(simYhats)
dim(simYhats)

predictionFrame <- data.frame(someScenarios,t(simYhats))

bulk <- apply(predictionFrame[1:101,11:1010],1,sort) 
bulk <- t(bulk) 
bulk2 <- bulk[,27:975]
bulk3 <- cbind(predictionFrame[,1:10],bulk2)

longFrame <- melt(predictionFrame, id.vars = 
	colnames(someScenarios))
head(longFrame)

test1 <- seq(-0.4,0.6,0.01)
test2 <- tapply(longFrame$value,longFrame$Var2,mean)
variable <- rep(1,101)
test <- data.frame(test1,test2,variable)

zp1 <- ggplot (data = longFrame,
		aes(x = Var2, y = value, group= paste(variable)))		
zp1 <- zp1 + geom_line(alpha=I(1/sqrt(nSims)))				
zp1 <- zp1 + geom_line(x=test1,y=test2,data=test,
		colour="white",size=1) # only for line at mean
zp1 <- zp1 + scale_x_continuous('Judicial independence',
		breaks=c(-0.4,-0.2,0,0.2,0.4,0.6),labels=c(
		"0.0","0.2","0.4","0.6","0.8","1.0"))
#zp1 <- zp1 + xlim(0,1)
zp1 <- zp1 + scale_y_continuous('',limits=c(0.2,1))
#zp1 <- zp1 + theme(plot.title="")
zp1 <- zp1 + theme_bw()
zp1 <- zp1 + theme(plot.title=element_text(size=25))
zp1 <- zp1 + theme(axis.title.y = element_text(size=24,angle=90),
          axis.text.y=element_text(size=20))
zp1 <- zp1 + theme(axis.title.x = element_text(size=24,vjust=0),
          axis.text.x  = element_text(size=20))
zp1  
        
##
# MODEL FOUR
##
res.model4 <- glm(ptfateD ~ SLJI + chga_demo, data3,family=binomial)

myModel <- res.model4
nSims <- 1000
someScenarios <- expand.grid(1,
	seq(-0.4,0.6,0.01),
	democ)
set.seed(1)	
simDraws <- mvrnorm(nSims,coef(myModel),vcov(myModel))
simYhats <- simDraws %*% t(someScenarios)
simYhats <- plogis(simYhats)
dim(simYhats)

predictionFrame <- data.frame(someScenarios,t(simYhats))

bulk <- apply(predictionFrame[1:101,4:1003],1,sort) 
bulk <- t(bulk) 
bulk2 <- bulk[,27:975]
bulk3 <- cbind(predictionFrame[,1:3],bulk2)

longFrame <- melt(predictionFrame, id.vars = 
	colnames(someScenarios))
head(longFrame)

test1 <- seq(-0.4,0.6,0.01)
test2 <- tapply(longFrame$value,longFrame$Var2,mean)
variable <- rep(1,101)
test <- data.frame(test1,test2,variable)

zp1 <- ggplot (data = longFrame,
		aes(x = Var2, y = value, group= paste(variable)))		
zp1 <- zp1 + geom_line(alpha=I(1/sqrt(nSims)))				
zp1 <- zp1 + geom_line(x=test1,y=test2,data=test,
		colour="white",size=1) # only for line at mean
zp1 <- zp1 + scale_x_continuous('',
		breaks=c(-0.4,-0.2,0,0.2,0.4,0.6),labels=c(
		"0.0","0.2","0.4","0.6","0.8","1.0"))
#zp1 <- zp1 + xlim(0,1)
zp1 <- zp1 + scale_y_continuous('Pr(Unpunished)',limits=c(0.2,1))
#zp1 <- zp1 + theme(plot.title="")
zp1 <- zp1 + theme_bw()
zp1 <- zp1 + theme(plot.title=element_text(size=25))
zp1 <- zp1 + theme(axis.title.y = element_text(size=24,angle=90),
          axis.text.y=element_text(size=20))
zp1 <- zp1 + theme(axis.title.x = element_text(size=24,vjust=0),
          axis.text.x  = element_text(size=20))
zp1          


##
## Descriptive stats
##

# Figure 1:
ggplot(data=fates, aes(x=ptfate2)) +  geom_histogram(binwidth=1, colour="black", fill="white")  +
		scale_y_continuous('Observations') + theme(legend.position=c(1,1)) +
		theme(legend.justification=c(1,1)) + theme(legend.title=element_text(size=18)) +
		theme(legend.text=element_text(size=14)) + 
		scale_x_discrete('',limits=c("Unpunished","Exiled","Imprisoned","Killed")) + 
		theme(axis.title.x = element_text(size=18),axis.text.x  = element_text(size=14)) +
		theme(axis.title.y = element_text(size=18,angle=90),axis.text.y=element_text(size=14))

# Figure 2:
form1 <- cname ~ year.x + ptfate2 + ptfateD + SLJI + LJI
data1 <- extractdata(form1,fates,na.rm=T)
plot <- ggplot(data1, aes(x=LJI)) + 
		geom_histogram(binwidth=0.02, colour="black", fill="white") +
		scale_y_continuous('Observations') + 
		scale_x_continuous('Judicial independence') +
		theme(axis.title.x = element_text(size=18),
          axis.text.x  = element_text(size=14)) +
        theme(axis.title.y = element_text(size=18,angle=90),
          axis.text.y=element_text(size=14)) 
plot          


##
## Appendix replication
##

##
# ROC PLOT
##
library(ROCR)
library(verification)
library(pROC)

roc1 <- roc(data1$ptfateD,res.model1$fitted)
#Area under the curve: 0.83
roc2 <- roc(data2$ptfateD,res.model2$fitted)
#roc1 <- roc under the curve: 0.85
roc3 <- roc(data3$ptfateD,res.model3$fitted)
#Area under the curve: 0.867
roc4 <- roc(data3$ptfateD,res.model4$fitted)
#Area under the curve: 0.850

pest <- c(roc1$auc,roc2$auc,roc3$auc,roc4$auc)
# The following is an alternate way, as the ci.auc command gives lower bounds [1],
# upper bounds [3], and *point estimates* [2]
#pest <- c(ci.auc(roc1)[2],ci.auc(roc2)[2],ci.auc(roc3)[2],ci.auc(roc4)[2],
#		ci.auc(roc5)[2],ci.auc(roc6)[2],ci.auc(roc7)[2],ci.auc(roc8)[2],ci.auc(roc9)[2])
		
lest <- c(ci.auc(roc1)[1],ci.auc(roc2)[1],ci.auc(roc3)[1],ci.auc(roc4)[1])
uest <- c(ci.auc(roc1)[3],ci.auc(roc2)[3],ci.auc(roc3)[3],ci.auc(roc4)[3])

# Below flips it so models are listed on the Y axis. 
# Note the need to reverse the order and re-label the y-axis,
# as normally it wants to put 1 directly above the x-axis

dat1 <- data.frame(x=seq(4,1,-1),pest=rev(pest),lower=rev(lest),upper=rev(uest))
plot <- ggplot(data=dat1, aes(x=x)) + 
		geom_segment(aes(y=seq(1,4,1),x=lower,yend=seq(1,4,1),xend=upper),size=1)
plot <- plot + geom_point(aes(x=pest,y=seq(1,4,1)),shape=19,size=4)
plot <- plot + scale_x_continuous('Area under the ROC curve',limits=c(0.75,0.95)) 
 # theme(legend.position=c(1,1)) +
# theme(legend.justification=c(1,1)) +
# theme(legend.title=element_text(size=18)) +
# theme(legend.text=element_text(size=14)) +
#theme(legend.background = theme_rect(fill="gray95")) +
plot <- plot + theme_bw()
plot <- plot + scale_y_continuous('',breaks=c(1,2,3,4),
		labels=c("Model 4","Model 3","Model 2","Model 1"))
plot <- plot + theme(axis.title.x = element_text(size=18),
           axis.text.x  = element_text(size=16))
plot <- plot + theme(axis.title.y = element_text(size=18),
           axis.text.y  = element_text(size=16))
plot


##
# CROSS VALIDATION
##
# Model 1:
set.seed(3)
infitset <- sample(rownames(data1),size=dim(data1)[[1]]/2) 
#assign rows to a fit set (fit being the set we "train" our model on)
totalset <- rownames(data1) #make a total set
intestset <- setdiff(totalset,infitset) #differentiate
fitset <- data1[infitset,] #make the fit set
testset <- data1[intestset,] #make the test set (to test the coef on)
fit1 <- glm(ptfateD ~ SLJI,data = fitset,family="binomial")

preds <- predict.glm(fit1,testset,type="response")
roc.area(data1$ptfateD,res.model1$fitted)[1] # 0.83
roc.area(fitset$ptfateD,fit1$fitted)[1] # 0.83
roc.area(testset$ptfateD,preds)[1]  # 0.83

# Model 2:
set.seed(3)
infitset <- sample(rownames(data2),size=dim(data2)[[1]]/2) 
#assign rows to a fit set (fit being the set we "train" our model on)
totalset <- rownames(data2) #make a total set
intestset <- setdiff(totalset,infitset) #differentiate
fitset <- data2[infitset,] #make the fit set
testset <- data2[intestset,] #make the test set (to test the coef on)
fit2 <- glm(ptfateD ~ SLJI + chga_demo + Slog_gdp,data = fitset,family="binomial")

preds <- predict.glm(fit2,testset,type="response")
roc.area(fitset$ptfateD,fit2$fitted)[1] # 0.86 
roc.area(testset$ptfateD,preds)[1]  # 0.85
roc.area(data2$ptfateD,res.model2$fitted)[1] # 0.85

# Model 3:
set.seed(3)
infitset <- sample(rownames(data3),size=dim(data3)[[1]]/2) 
#assign rows to a fit set (fit being the set we "train" our model on)
totalset <- rownames(data3) #make a total set
intestset <- setdiff(totalset,infitset) #differentiate
fitset <- data3[infitset,] #make the fit set
testset <- data3[intestset,] #make the test set (to test the coef on)
fit3 <- glm(ptfateD ~ SLJI + chga_demo + Slog_gdp + major_war + 
	Syears_democ + as.factor(entry) + dpi_system,data = fitset,family="binomial")

preds <- predict.glm(fit3,testset,type="response")
roc.area(fitset$ptfateD,fit3$fitted)[1] # 0.87
roc.area(testset$ptfateD,preds)[1]  # 0.86
roc.area(data3$ptfateD,res.model3$fitted)[1] # 0.87

# Model 4:
set.seed(3)
infitset <- sample(rownames(data3),size=dim(data3)[[1]]/2) 
#assign rows to a fit set (fit being the set we "train" our model on)
totalset <- rownames(data3) #make a total set
intestset <- setdiff(totalset,infitset) #differentiate
fitset <- data3[infitset,] #make the fit set
testset <- data3[intestset,] #make the test set (to test the coef on)
fit4 <- glm(ptfateD ~ SLJI + chga_demo,data = fitset,family="binomial")

preds <- predict.glm(fit4,testset,type="response")
roc.area(fitset$ptfateD,fit4$fitted)[1] # 0.85
roc.area(testset$ptfateD,preds)[1]  # 0.85
roc.area(data3$ptfateD,res.model4$fitted)[1] # 0.85

##
# Excluding the killed category
##

fates1$ptfateDNK <- fates1$ptfateD
fates1$ptfateDNK[fates1$ptfate2==4] <- NA

form1 <- cname ~ year.x + ptfate2 + ptfateDNK + SLJI + leader
form2 <- cname ~ year.x + ptfate2 + ptfateDNK + ptfateD + SLJI + chga_demo + Slog_gdp
form3 <- cname ~ year.x + ptfate2 + ptfateDNK + ptfateD + SLJI + chga_demo + 
		Slog_gdp + dpi_system + major_war + entry + Syears_democ

data1 <- extractdata(form1,fates1,na.rm=T)
data2 <- extractdata(form2,fates1,na.rm=T)	
data3 <- extractdata(form3,fates1,na.rm=T)	

out1k <- glm(ptfateDNK ~ SLJI,data1,family=binomial)
out2k <- glm(ptfateDNK ~ SLJI + chga_demo + Slog_gdp,data2,family=binomial)
out3k <- glm(ptfateDNK ~ SLJI + chga_demo + Slog_gdp  + major_war +
			Syears_democ + as.factor(entry) + dpi_system,
			data3,family=binomial)		
out4k <- glm(ptfateDNK ~ SLJI + chga_demo,data3,family=binomial)
m1 <- mtable(out1k,out2k,out3k,out4k,			
		summary.stats=c("BIC","N"))			
m2 <- (m1 <- relabel(m1,
		"Ranef - cname" = "Random effect: country",
		"SLJI" = "Judicial independence",
		"chga_demo: 1. Democracy/0. Dictatorship" = "Democracy",
		"Slog_gdp" = "log(GDP/capita)",
		"major_war" = "Conflict",
		"Swdi_gdpgr" = "Economic growth",
		"Syears_democ" = "Years democratic",
		"as.factor(entry): 1/0" = "Entry by irregular means",
		"as.factor(entry): 2/0" = "Entry by foreign imposition",
		"dpi_system: 1. Strong president elected by assembly/0. Direct presidential" =
		"President elected by assembly",
		"dpi_system: 2. Parliamentary/0. Direct presidential" = "Parliamentary system"
		))
m2

##
# Considering who replaces
##


opps <- read.csv("replication_data_2.csv")
form1o <- cname ~ year.x + ptfateD + SLJI + opposition
form2o <- cname ~ year.x + ptfateD + SLJI + chga_demo + Slog_gdp + opposition
form3o <- form1o <- cname ~ year.x + ptfateD + SLJI + chga_demo + Slog_gdp +
				major_war + Syears_democ +
				as.factor(entry) + dpi_system + opposition
data1o <- extractdata(form1o,opps,na.rm=TRUE) 
data2o <- extractdata(form2o,opps,na.rm=TRUE) 
data3o <- extractdata(form3o,opps,na.rm=TRUE) 

## when opposition=0 same party replaces, when =1 opponents replace

out1o <- glm(ptfateD ~ SLJI ,
				data=data1o,subset=opposition==0,family=binomial)
out1no <- glm(ptfateD ~ SLJI ,
				data=data1o,subset=opposition==1,family=binomial)

out2o <- glm(ptfateD ~ SLJI + chga_demo + Slog_gdp,
				data=data2o,subset=opposition==0,family=binomial)
out2no <- glm(ptfateD ~ SLJI + chga_demo + Slog_gdp,
				data=data2o,subset=opposition==1,family=binomial)
				
out3o <- glm(ptfateD ~ SLJI + chga_demo + Slog_gdp + major_war + Syears_democ +
				as.factor(entry) + dpi_system,
				data=data3o,subset=opposition==0,family=binomial)
out3no <- glm(ptfateD ~ SLJI + chga_demo + Slog_gdp + major_war + Syears_democ +
				as.factor(entry) + dpi_system,
				data=data3o,subset=opposition==1,family=binomial)

out4o <- glm(ptfateD ~ SLJI + chga_demo,
				data=data3o,subset=opposition==0,family=binomial)
out4no <- glm(ptfateD ~ SLJI + chga_demo,
				data=data3o,subset=opposition==1,family=binomial)
								
								
m1 <- mtable(out1o,out1no,out2o,out2no,out3o,out3no,out4o,out4no,
		summary.stats=c("BIC","N"))			
m2 <- (m1 <- relabel(m1,
		"Ranef - cname" = "Random effect: country",
		"SLJI" = "Judicial independence",
		"chga_demo: 1. Democracy/0. Dictatorship" = "Democracy",
		"Slog_gdp" = "log(GDP/capita)",
		"major_war" = "Conflict",
		"Swdi_gdpgr" = "Economic growth",
		"Syears_democ" = "Years democratic",
		"as.factor(entry): 1/0" = "Entry by irregular means",
		"as.factor(entry): 2/0" = "Entry by foreign imposition",
		"dpi_system: 1. Strong president elected by assembly/0. Direct presidential" =
		"President elected by assembly",
		"dpi_system: 2. Parliamentary/0. Direct presidential" = "Parliamentary system"
		))
m2	
			

##
# HIERARCHICAL LOGIT
##
form1 <- cname ~ year.x + ptfate2 + ptfateD + SLJI 
data1 <- extractdata(form1,fates1,na.rm=T)
form2 <- cname ~ year.x + ptfate2 + ptfateD + SLJI + chga_demo + Slog_gdp			
data2 <- extractdata(form2,fates1,na.rm=T)
form3 <- cname ~ year.x + ptfate2 + ptfateD + SLJI + chga_demo + Slog_gdp +
 	major_war + Syears_democ + entry + dpi_system	
data3 <- extractdata(form3,fates1,na.rm=T) #model 4 uses the same data as model 3

library(lme4)
h.out1 <- glmer(ptfateD ~ SLJI + (1|cname),data1,family="binomial",REML=F)
h.out2 <- glmer(ptfateD ~ SLJI + chga_demo + Slog_gdp + 
			(1|cname),data2,family="binomial",REML=F)
h.out3 <- glmer(ptfateD ~ SLJI + chga_demo + Slog_gdp  + major_war +
			Syears_democ + as.factor(entry) + dpi_system +
			(1|cname),data3,family="binomial",REML=F)
h.out4 <- glmer(ptfateD ~ SLJI + chga_demo + (1|cname),data3,family="binomial",REML=F)
		
m1 <- mtable(h.out1,h.out2,h.out3,h.out4,
		summary.stats=c("BIC","N"))			
m2 <- (m1 <- relabel(m1,
		"Ranef - cname" = "Random effect: country",
		"SLJI" = "Judicial independence",
		"chga_demo: 1. Democracy/0. Dictatorship" = "Democracy",
		"Slog_gdp" = "log(GDP/capita)",
		"major_war" = "Conflict",
		"Swdi_gdpgr" = "Economic growth",
		"Syears_democ" = "Years democratic",
		"as.factor(entry): 1/0" = "Entry by irregular means",
		"as.factor(entry): 2/0" = "Entry by foreign imposition",
		"dpi_system: 1. Strong president elected by assembly/0. Direct presidential" =
		"President elected by assembly",
		"dpi_system: 2. Parliamentary/0. Direct presidential" = "Parliamentary system"
		))
m2		
			
##
# ORDERED LOGIT
##
fates1$ptfate3 <- 0
fates1$ptfate3[fates1$ptfate2==4] <- 1
fates1$ptfate3[fates1$ptfate2==3] <- 2
fates1$ptfate3[fates1$ptfate2==2] <- 3
fates1$ptfate3[fates1$ptfate2==1] <- 4

form1 <- cname ~ year.x + ptfate3 + ptfateD + SLJI
form2 <- cname ~ year.x + ptfate3 + ptfateD + SLJI + chga_demo + Slog_gdp			
form3 <- cname ~ year.x + ptfate3 + ptfateD + SLJI + chga_demo + Slog_gdp + 
		dpi_system + major_war + entry + Syears_democ
		
data1 <- extractdata(form1,fates1,na.rm=T)
data2 <- extractdata(form2,fates1,na.rm=T)
data3 <- extractdata(form3,fates1,na.rm=T)	



p.out1 <- polr(as.factor(ptfate3) ~ SLJI,data1)
p.out2 <- polr(as.factor(ptfate3) ~ SLJI + chga_demo + Slog_gdp,data2)
p.out3 <- polr(as.factor(ptfate3) ~ SLJI + chga_demo + Slog_gdp +
	major_war + Syears_democ + as.factor(entry) + dpi_system,data3)
p.out4 <- polr(as.factor(ptfate3) ~ SLJI + chga_demo,data3)
	

m1 <- mtable(p.out1,p.out2,p.out3,p.out4,			
		summary.stats=c("BIC","N"))			
m2 <- (m1 <- relabel(m1,
		"Ranef - cname" = "Random effect: country",
		"SLJI" = "Judicial independence",
		"chga_demo: 1. Democracy/0. Dictatorship" = "Democracy",
		"Slog_gdp" = "log(GDP/capita)",
		"major_war" = "Conflict",
		"Swdi_gdpgr" = "Economic growth",
		"Syears_democ" = "Years democratic",
		"as.factor(entry): 1/0" = "Entry by irregular means",
		"as.factor(entry): 2/0" = "Entry by foreign imposition",
		"dpi_system: 1. Strong president elected by assembly/0. Direct presidential" =
		"President elected by assembly",
		"dpi_system: 2. Parliamentary/0. Direct presidential" = "Parliamentary system"
		))
m2		
		
## 
# ORDERED LOGIT FIGURES
##		

# Model 1:
beta2 <- p.out1$coef
tau2 <- p.out1$zeta
AB <- seq(-0.4,0.6,0.01)
X <- cbind(AB)

logit.prob <- function(eta){exp(eta)/(1+exp(eta))}
pA <- logit.prob(tau2[1] - X %*% beta2)
pB <- logit.prob(tau2[2] - X %*% beta2) - logit.prob(tau2[1] - X %*% beta2)
pC <- logit.prob(tau2[3] - X %*% beta2) - logit.prob(tau2[2] - X %*% beta2)
pD <- 1.0 - logit.prob(tau2[3] - X %*% beta2)

hmm <- cbind(pA,pB,pC,pD)
dat <- data.frame(x=seq(-0.4,0.6,0.01),killed=hmm[,1],
			imprisoned=hmm[,2],
			exiled=hmm[,3],
			unpunished=hmm[,4]
			)

plot <- ggplot(dat, aes(x=x, y= Probability))
plot <- plot + geom_line(aes(y=killed),linetype=1,size=1,colour="gray70")
plot <- plot + geom_line(aes(y=imprisoned),linetype=3,size=1.25)
plot <- plot + geom_line(aes(y=exiled),linetype=2,size=1,colour="gray70")
plot <- plot + geom_line(aes(y=unpunished),linetype=1,size=1)
plot <- plot + scale_x_continuous('',
		breaks=c(-0.4,-0.2,0,0.2,0.4,0.6),labels=c(
		"0.0","0.2","0.4","0.6","0.8","1.0"))
plot <- plot + scale_y_continuous('E(Unpunished)',limits=c(0,1))
plot <- plot + theme_bw()
plot <- plot + theme(plot.title=element_text(size=25))
plot <- plot + theme(axis.title.y = element_text(size=24,angle=90),
          axis.text.y=element_text(size=20))
plot <- plot + theme(axis.title.x = element_text(size=18,vjust=0),
          axis.text.x  = element_text(size=20))
plot

# Model 2:
p.out2 <- polr(as.factor(ptfate3) ~ SLJI + chga_demo + Slog_gdp, 
    data = data2)
beta2 <- p.out2$coef
tau2 <- p.out2$zeta
AB <- seq(-0.4,0.6,0.01)
X <- cbind(AB,
		   0.6,
		   mean(data3$Slog_gdp)) # 
		   
logit.prob <- function(eta){exp(eta)/(1+exp(eta))}
pA <- logit.prob(tau2[1] - X %*% beta2)
pB <- logit.prob(tau2[2] - X %*% beta2) - logit.prob(tau2[1] - X %*% beta2)
pC <- logit.prob(tau2[3] - X %*% beta2) - logit.prob(tau2[2] - X %*% beta2)
pD <- 1.0 - logit.prob(tau2[3] - X %*% beta2)

hmm <- cbind(pA,pB,pC,pD)
dat <- data.frame(x=seq(-0.4,0.6,0.01),killed=hmm[,1],
			imprisoned=hmm[,2],
			exiled=hmm[,3],
			unpunished=hmm[,4]
			)

plot <- ggplot(dat, aes(x=x, y= Probability))
plot <- plot + geom_line(aes(y=killed),linetype=1,size=1,colour="gray70")
plot <- plot + geom_line(aes(y=imprisoned),linetype=3,size=1.25)
plot <- plot + geom_line(aes(y=exiled),linetype=2,size=1,colour="gray70")
plot <- plot + geom_line(aes(y=unpunished),linetype=1,size=1)
plot <- plot + scale_x_continuous('Judicial independence',
		breaks=c(-0.4,-0.2,0,0.2,0.4,0.6),labels=c(
		"0.0","0.2","0.4","0.6","0.8","1.0"))
plot <- plot + scale_y_continuous('',limits=c(0,1))
plot <- plot + theme_bw()
plot <- plot + theme(plot.title=element_text(size=25))
plot <- plot + theme(axis.title.y = element_text(size=24,angle=90),
          axis.text.y=element_text(size=20))
plot <- plot + theme(axis.title.x = element_text(size=24,vjust=0),
          axis.text.x  = element_text(size=20))
plot

# Model 3:

beta2 <- p.out3$coef
tau2 <- p.out3$zeta
AB <- seq(-0.4,0.6,0.01)
X <- cbind(AB,
		   0.6,
		   mean(data8$Slog_gdp),
		   mean(data8$major_war),
		   mean(data8$Syears_democ),
		   0,
		   0,
		   0,
		   0)
		   
logit.prob <- function(eta){exp(eta)/(1+exp(eta))}
pA <- logit.prob(tau2[1] - X %*% beta2)
pB <- logit.prob(tau2[2] - X %*% beta2) - logit.prob(tau2[1] - X %*% beta2)
pC <- logit.prob(tau2[3] - X %*% beta2) - logit.prob(tau2[2] - X %*% beta2)
pD <- 1.0 - logit.prob(tau2[3] - X %*% beta2)

hmm <- cbind(pA,pB,pC,pD)
dat <- data.frame(x=seq(-0.4,0.6,0.01),killed=hmm[,1],
			imprisoned=hmm[,2],
			exiled=hmm[,3],
			unpunished=hmm[,4]
			)

plot <- ggplot(dat, aes(x=x, y= Probability))
plot <- plot + geom_line(aes(y=killed),linetype=1,size=1,colour="gray70")
plot <- plot + geom_line(aes(y=imprisoned),linetype=3,size=1.25)
plot <- plot + geom_line(aes(y=exiled),linetype=2,size=1,colour="gray70")
plot <- plot + geom_line(aes(y=unpunished),linetype=1,size=1)
plot <- plot + scale_x_continuous('Judicial independence',
		breaks=c(-0.4,-0.2,0,0.2,0.4,0.6),labels=c(
		"0.0","0.2","0.4","0.6","0.8","1.0"))
plot <- plot + scale_y_continuous('',limits=c(0,1))
plot <- plot + theme_bw()
plot <- plot + theme(plot.title=element_text(size=25))
plot <- plot + theme(axis.title.y = element_text(size=24,angle=90),
          axis.text.y=element_text(size=20))
plot <- plot + theme(axis.title.x = element_text(size=24,vjust=0),
          axis.text.x  = element_text(size=20))
plot
			
# Model 4:

beta2 <- p.out4$coef
tau2 <- p.out4$zeta
AB <- seq(-0.4,0.6,0.01)
X <- cbind(AB,
		   0.6)
		   
logit.prob <- function(eta){exp(eta)/(1+exp(eta))}
pA <- logit.prob(tau2[1] - X %*% beta2)
pB <- logit.prob(tau2[2] - X %*% beta2) - logit.prob(tau2[1] - X %*% beta2)
pC <- logit.prob(tau2[3] - X %*% beta2) - logit.prob(tau2[2] - X %*% beta2)
pD <- 1.0 - logit.prob(tau2[3] - X %*% beta2)

hmm <- cbind(pA,pB,pC,pD)
dat <- data.frame(x=seq(-0.4,0.6,0.01),killed=hmm[,1],
			imprisoned=hmm[,2],
			exiled=hmm[,3],
			unpunished=hmm[,4]
			)

plot <- ggplot(dat, aes(x=x, y= Probability))
plot <- plot + geom_line(aes(y=killed),linetype=1,size=1,colour="gray70")
plot <- plot + geom_line(aes(y=imprisoned),linetype=3,size=1.25)
plot <- plot + geom_line(aes(y=exiled),linetype=2,size=1,colour="gray70")
plot <- plot + geom_line(aes(y=unpunished),linetype=1,size=1)
plot <- plot + scale_x_continuous('',
		breaks=c(-0.4,-0.2,0,0.2,0.4,0.6),labels=c(
		"0.0","0.2","0.4","0.6","0.8","1.0"))
plot <- plot + scale_y_continuous('E(Unpunished)',limits=c(0,1))
plot <- plot + theme_bw()
plot <- plot + theme(plot.title=element_text(size=25))
plot <- plot + theme(axis.title.y = element_text(size=24,angle=90),
          axis.text.y=element_text(size=20))
plot <- plot + theme(axis.title.x = element_text(size=18,vjust=0),
          axis.text.x  = element_text(size=20))
plot
		
##
# MULTINOMIAL LOGIT
##

form2 <- cname ~ year.x + ptfate2 + ptfateD + SLJI + 
		chga_demo + Slog_gdp + LJI	+ ptfate3 + democ.cont
data2 <- extractdata(form2,fates1,na.rm=T)

data2$ptfate3.cat <- data2$ptfate2
data2$ptfate3.cat[data2$ptfate3.cat==4] <- 3 # collapsing the killed and imprisoned categories for Model 2g

library(nnet)
mlogit.result <- multinom(as.factor(ptfate2) ~ SLJI + democ.cont + 
	Slog_gdp,Hess=TRUE,data=data2) # Model 2f
mlogit.3cat <- multinom(as.factor(ptfate3.cat) ~ SLJI + democ.cont + Slog_gdp,Hess=TRUE,data=data2)	# Model 2g

##
# MLOGIT FIGURES
##

# figure for Model 2f
predictors <- expand.grid(SLJI=seq(-0.4,0.6,0.01),
						democ.cont=mean(data2$democ.cont),
						Slog_gdp=mean(data2$Slog_gdp))
p.fit <- predict(mlogit.result,predictors,type='probs')
huh<- data.frame(predictors,p.fit)
p1 <- huh[,4]
p2 <- huh[,5]
p3 <- huh[,6]
p4 <- huh[,7]

dat1 <- data.frame(x=seq(-0.4,0.6,0.01),unpunished=p1,exiled=p2,imprisoned=p3,killed=p4)

plot <- ggplot(dat1, aes(x=x, y= Probability))
plot <- plot + geom_line(aes(y=killed),linetype=1,size=1,colour="gray70")
plot <- plot + geom_line(aes(y=imprisoned),linetype=3,size=1.25)
plot <- plot + geom_line(aes(y=exiled),linetype=2,size=1,colour="gray70")
plot <- plot + geom_line(aes(y=unpunished),linetype=1,size=1)
plot <- plot + scale_x_continuous('Judicial independence',
		breaks=c(-0.4,-0.2,0,0.2,0.4,0.6),labels=c(
		"0.0","0.2","0.4","0.6","0.8","1.0"))
plot <- plot + scale_y_continuous('E(Unpunished)',limits=c(0,1))
plot <- plot + theme_bw()
plot <- plot + theme(plot.title=element_text(size=25))
plot <- plot + theme(axis.title.y = element_text(size=24,angle=90),
          axis.text.y=element_text(size=20))
plot <- plot + theme(axis.title.x = element_text(size=24,vjust=0),
          axis.text.x  = element_text(size=20))
plot

# Figure for Model 2g

predictors <- expand.grid(SLJI=seq(-0.4,0.6,0.01),
						democ.cont=mean(data2$democ.cont),
						Slog_gdp=mean(data2$Slog_gdp))
p.fit <- predict(mlogit.3cat,predictors,type='probs')
huh<- data.frame(predictors,p.fit)
p1 <- huh[,4]
p2 <- huh[,5]
p3 <- huh[,6]

dat1 <- data.frame(x=seq(-0.4,0.6,0.01),unpunished=p1,exiled=p2,imprisonedkilled=p3)

plot <- ggplot(dat1, aes(x=x, y= Probability))
plot <- plot + geom_line(aes(y=imprisonedkilled),linetype=3,size=1.25)
plot <- plot + geom_line(aes(y=exiled),linetype=2,size=1,colour="gray70")
plot <- plot + geom_line(aes(y=unpunished),linetype=1,size=1)
plot <- plot + scale_x_continuous('Judicial independence',
		breaks=c(-0.4,-0.2,0,0.2,0.4,0.6),labels=c(
		"0.0","0.2","0.4","0.6","0.8","1.0"))
plot <- plot + scale_y_continuous('',limits=c(0,1))
plot <- plot + theme_bw()
plot <- plot + theme(plot.title=element_text(size=25))
plot <- plot + theme(axis.title.y = element_text(size=24,angle=90),
          axis.text.y=element_text(size=20))
plot <- plot + theme(axis.title.x = element_text(size=24,vjust=0),
          axis.text.x  = element_text(size=20))
plot



		